
program define xtprval , rclass
version 10.0
qui {
syntax varlist [,  SEED(int 0)] AT(numlist max=2) REST(string) Level(cilevel)

local ecmd = e(cmdline)
gettoken modname ecmd : ecmd
tempvar touse
gen `touse' = e(sample)

/*
if "`modname'"~="xtlogit" & "`modname'"~="xtprobit" {
di  as error "The command only works with xtlogit and xtprobit"
exit 198
}
*/
mat define bm=get(_b)
mat  v = e(V)
mat vm = vecdiag(cholesky(v))

local n: colname bm

local ck: list varlist & n
if "`ck'"~="`varlist'" {
nois di as error "Variable `varlist' was not in the model."
exit 198
}

local ct : word count `at'
// Differences for PA and RE models
// with RE need to set ln_sigma

scalar c = colsof(bm)
mat vals = J(1,colsof(bm),0)
mat coln vals = `n'

mat vals[1,colsof(bm)]=1
local ivnum = colsof(bm)-1
 forval i = 1/`ivnum' {
 local v : word `i' of `n'
 qui sum `v' if `touse'
 mat vals[1,`i'] = r(`rest')
 }
 if `ct'==1 {
 mat vals[1, colnumb(vals,"`varlist'")]=`at'
 }
 if `ct'==2 {
 local at1: word 1 of `at'
  mat vals[1, colnumb(vals,"`varlist'")]=`at1'
  mat vals2=vals
   local at2: word 2 of `at'
  mat vals2[1, colnumb(vals,"`varlist'")]=`at2'
}
 
 
* gen s = .
mata b = st_matrix("bm")
mata v = st_matrix("v")

mata bsims =b :+ (invnormal(uniform(1000,cols(v)))*cholesky(v)')

mata vals = st_matrix("vals")
mata s1 = normal(bsims* vals')

foreach 9 in _s1 _s2 _fd {
capture confirm new variable `9'
if _rc~=0 {
di "Dropping existing variable `9'"
drop `9'
}
}

mata ss1 = st_addvar("float","_s1")
/// getmata s1 , force replace
mata st_store ((1,rows(s1)),ss1,s1)
if `ct'==2 {
mata vals2 = st_matrix("vals2")
mata s2 = normal(bsims* vals2')
///getmata s2 , force replace
mata ss2 = st_addvar("float","_s2")
mata st_store((1,rows(s2)),ss2,s2)
capture drop fd
gen _fd = _s2-_s1
}
if `ct'==1 {
local s _s1
summ `s'
local mu =r(mean)
tempname ll ul
local `ll' = round((100-`level')/2,.01)
local `ul' = round((100+`level')/2,.01)
di "``ll''" " ``ul''"
_pctile `s' , p(``ll'' ``ul'')
 local `s'lo=r(r1)
local `s'hi=r(r2)
return scalar s1=`mu'
return scalar `s'lo=``s'lo'
return scalar `s'hi=``s'hi'
nois {
di as text "{hline 33}{c TT}{hline 52}"
di as text   "        Variable                 {c |}  Prediction        [`level'% Conf. Interval]"
di as text "{hline 33}{c +}{hline 52}"
di as result "        `varlist'   `at' "   _col(34)  "{c |}"           as result %9.0g     `s' "         "   as result %9.0g  ``s'lo'  "   "   as result %9.0g  ``s'hi'    "
di as text "{hline 33}{c BT}{hline 52}"
}

}
if `ct'==2 {
foreach s of varlist _s1 _s2 _fd {
qui summ `s'
local `s'mu =r(mean)
tempname ll ul
local `ll' = round((100-`level')/2,.01)
local `ul' = round((100+`level')/2,.01)
_pctile `s' , p(``ll'' ``ul'')
 local `s'lo=r(r1)
local `s'hi=r(r2)
return scalar `s'=``s'mu'
return scalar `s'lo=``s'lo'
return scalar `s'hi=``s'hi'
}
nois {
di as text "{hline 33}{c TT}{hline 52}"
di as text   "        Variable" _col(34) "{c |}  Prediction        [`level'% Conf. Interval]"
di as text "{hline 33}{c +}{hline 52}"
di as result "        `varlist'  `at1' "   _col(34)  "{c |}"           as result %9.0g     `_s1mu' "          "   as result %9.0g  `_s1lo'  "   "   as result %9.0g  `_s1hi'    "
di as result "        `varlist'  `at2' "   _col(34)  "{c |}"           as result %9.0g     `_s2mu' "          "   as result %9.0g  `_s2lo'  "   "   as result %9.0g  `_s2hi'    "
di as result "        Difference "         _col(34)  "{c |}"           as result %9.0g     `_fdmu' "          "   as result %9.0g  `_fdlo'  "   "   as result %9.0g  `_fdhi'    "

di as text "{hline 33}{c BT}{hline 52}"
}


}
}

end
